Statistical mechanics of ecosystem assembly 
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We introduce a toy model of ecosystem assembly for which we are able to map out all assembly pathways 
generated by external invasions. The model allows to display the whole phase space in the form of an assembly 
graph whose nodes are communities of species and whose directed links are transitions between them induced 
by invasions. We characterize the process as a finite Markov chain and prove that it exhibits a unique set of 
recurrent states (the endstate of the process), which is therefore resistant to invasions. This also shows that the 
endstate is independent on the assembly history. The model shares all features with standard assembly models 
reported in the literature, with the advantage that all observables can be computed in an exact manner 

PACS numbers: 87.23.Cc, 87.23.Kg, 89.75.Fb, 64.60.De 



Understanding why ecosystems are both stable and com- 
plex is still one of the open questions in Ecology [1]. At 
present, it is widely accepted that the dynamic process by 
which ecological communities are assembled is key to solve 
this puzzle. It also has profound implications for conser- 
vation — e.g. it can shed light on how biodiversity may re- 
cover after major crisis, or the influence it has on commu- 
nity stability. Although ecosystem assembly has been stud- 
ied experimentally 1 2, 3.1 , the bulk of studies are theoreti- 
cal ||i S S 0, II laKl. These assembly models are but 
idealizations of the complex processes taking place in real 
ecosystem assembly, which nevertheless implement the same 
mechanisms acting in the formation of real ecosystems 111 ill . 
In this respect these models are close in spirit to the gen- 
eral approach of statistical physics of devising oversimplified 
paradigms which provide insights into the real phenomena. 

Previous work on ecosystem assembly made stochastic re- 
alizations of sequential invasions based on a finite set of pos- 
sible invaders (known as "regional species pool" [7]) whose 
trophic interactions are determined in advance. These stud- 
ies conclude: (i) at the end of the process, a final endstate 
resistant to invasions is reached, which can be either a sin- 
gle community or a set involving more than one community 
between which the system fluctuates [8]; (ii) the average re- 
sistance of an ecosystem to be colonized increases in time, 
and (iii) species richness also increases in time. Then the as- 
sembly process favors increasing species richness as well as 
stability, understood as resistance to invasions. 

Successful as these models may be to provide insights into 
the formation of ecosystems, we still lack a global picture of 
the process. The reason is at least twofold. First, these model 
are defined as very complex stochastic processes not amenable 
to analytical treatment, so that one can only hope to simulate 
a limited set of realizations of the process and take averages 
on them. This leaves some questions open such as, for in- 
stance, whether the endstate depends on the assembly history. 
Althoiigh most evidence points to the uniqueness of this end- 
state illll . it is still a matter of discussion j 1 2tl . Secondly, the 
size of the typical species pool is always small, so that it has 



been argued III 311 that the exhaustion of good invaders might 
justify the increasing resistance to invasion. 

Our aim in this Letter is to propose a toy version of an as- 
sembly model which, despite the styhzed communities it deals 
with, still contains all the ingredients of standard assembly 
models (something like the 'Tsing model" of ecosystem as- 
sembly) and lacks some of its limitations (Uke the finiteness 
of the regional species pool). Its major advantage over them 
is that it allows to map out the whole "phase space" of the 
system, as well as to analyze all assembly pathways, some- 
thing that permits us to obtain exact results, to explore the 
set of parameters characterizing different regimes, and to an- 
swer questions out of reach of standard assembly models. The 
model recovers all results mentioned above and provides new 
ones — the independence of the endstate from the assembly 
history being the most prominent among them. 

In what follows ecosystem will refer to the system as a 
whole, whereas community will stand for any particular col- 
lection of species, i.e. any realization of the ecosystem. 

Two are the ingredients that assembly models build upon: A 
deterministic population dynamics for the species forming the 
community, and a stochastic invasion process. For the former 
we essentially adopt a mean-field-like Lotka-Volterra model 
recently proposed to study coexistence in competing commu- 
nities and trophic level organization in food webs j 14, 15i[l6ll 
(these kind of models are called neutral in Ecology 11711 . al- 
though neutral models assume a stochastic population dynam- 
ics, unlike the one defined here). In this model, communities 
are arranged in trophic levels, and species at level £ are as- 
sumed to feed only on species at level £ — \, and on all of 
them alike (mean field). Although omnivory (i.e. feeding from 
lower levels) can be accounted for, it is still a matter of debate 
how common it is |,18J. so we disregard it. Accordingly, if nf 
denotes the population of species i at level I < i < L, 
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where = ""A^ being the number of species at level 

i, 7+ controls the amount of energy available to reproduction 
for each predation event, and 7_ takes into account the mean 
damage caused by predation to prey reproduction. Following 
the standard rule of efficiency on upwards energy transport 
to the next trophic level ifioll . we assume that these two con- 
stants are related by 7+ = 0.l7_. Direct interspecific com- 
petition is measured by p < 1, while intraspecific competi- 
tion is set to one to fix the population scale [16]. We regard 
all these species as consumers, and so they have a death rate 
a > 0. For simplicity all parameters are assumed to be the 
same for all species (this assumption is harmless because the 
system has been shown to be structurally stable against vari- 
ations of the constants LiSJ). The second of Eqs. ([T]) states 
that all species at the first level (basal species) predate on a 
single resource whose "amount" is given by uq. In the ab- 
sence of basal species this resource reaches a steady value of 
no ~ R, thus this constant determines the maximum amount 
of resource available to consumers. Finally, as real popu- 
lations cannot be arbitrarily small, it is required that extant 
species have a population > ric, the extinction threshold. 
If a population falls below this value it is considered extinct 
and removed from the community (see below). Low popula- 
tions are vulnerable against external variations or adverse mu- 
tations (20], and this stochastic effect is somehow accounted 
for in the deterministic dynamics ([TJ by the introduction of 
this viability condition. Calculations have been carried out 
with a = 1, nc = 1, 7- = 5, p = 0.3 and < i? < 1700. 
We have checked that variations of these parameters do not 
change the qualitative picture. 

Equations ([T]i have several equilibria, the main one being 
the interior equilibrium. In it all species of level I have the 
same population n}. These are the solution to the system 



-a =7-Af;+i + [1 + [si - 1)p\n;/si - j+N^_„ 
i?=iV*+7-iV*, 



(2) 



where £ — 1, . . . , L, Ng ~ sgn*^, and we have the constraints 
So = 1, sl+1 = 0. The remaining equilibria are obtained 
by setting to zero any subset of the populations and solving 
the linear system resulting from eliminating those variables. 
Therefore one only needs the solutions of the linear systems 
(|2|i for all possible choices of Si in order to characterize the 
dynamical equilibrium of this model. Interior equilibria are 
globally stable if, and only if, > for all £ = 0, ... , L, 
because 
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is a Lyapunov function 1I21I1 for ([T]l- Therefore, given the set 
of species numbers {si}f^i, the corresponding community is 
viable and stable if, and only if, rt^ > ric for all£ — 1, . . . ,L. 
Thus by solving Q for all choices of species numbers we can 
determine all viable and stable communities. Let denote 
the set of these communities. Although the total size of a 




FIG. 1: (Color online.) Assembly graph obtained for 7? — 140. It 
displays 130 communities with up to three trophic levels. The num- 
bers indicate how many species are in each level. Black arrows rep- 
resent accepted invasions; red arrows represent transitions inducing a 
species loss. The only node in blue (corresponding to the community 
{16, 4, 1}) represents the endstate of the assembly process. 



trophic level is not explicitly limited in any way, the extinction 
threshold ric imposes a constraint and therefore T is finite. 

Equilibrium communities will undergo invasions which 
change their composition. We shall assume that they take 
place at a uniform rate Two hypotheses underlying this 
model, as well as most assembly models, are: (i) the typical 
dynamical time is much smaller than so that communi- 
ties are never invaded during transients, and (ii) the popula- 
tion of invaders is small 12211 . Obviously, if the invasion rate 
is too high 1231 12411 or if invaders arrive with high popula- 
tions 1251], the assembly dynamics — and consequently the fi- 



nal endstates — will be seriously altered. These situations are 
outside the scope of this model. 

The invasion process goes as follows. Consider a commu- 
nity E E T, with L trophic levels, at its equilibrium point. 
Invasions can occur at any level t = 1, . . . , L + 1, so we ran- 
domly select i and add a new species, with initial population 
Tie, at level I of the community E. Because of the global sta- 
bility of our model, the extended community evolves to the 
equilibrium given by (|2]i with sg replaced by + 1. If this 
equilibrium is viable, then we will have a new community E' 
consisting of E plus the invader, and a transition will have oc- 
curred from E to E' . If, on the contrary, the population of the 
species of at least one level falls below ric, then extinctions 
will occur. It is unrealistic to extinguish the whole level, even 
though all its species have the same population. In reality, if 
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FIG. 2: (Color online.) Endstate of the assembly process obtained 
for R = 430. It contains 68 communities (the assembly graph con- 
tains 3060 communities). Black arrows represent accepted invasions, 
whereas yellow ones represent transitions with species loss (width is 
proportional to the number of extinctions). The diameter of the nodes 
is proportional to the stationary probability of its community. Labels 
provide the number of species at each level. 



many species are threatened, by chance one of them will be 
the first to become extinct. This fact may help the remaining 
species to survive. Accordingly we shall extinguish species in 
an inviable level as follows. As we are monitoring the whole 
trayectory of the system, we can detect the the moment when 
the first trophic level crossed Uc- At that point we remove 
one species from that level and restart the evolution from that 
point. We keep on removing one species at a time and restart- 
ing until the new resulting equilibrium becomes viable. Two 
things can thus happen: Either the first level to fall below ric is 
the invaded level, in which case the invader is simply rejected 
and no transition occurs, or it is another level that falls below 
Uc- In this case we end up with a new community E", and a 
transition will have occurred from E to E" . 

Starting from the empty community, 0, and considering 
for every community all possible invasions, we construct a 
directed graph, the assembly graph, Q, whose nodes are ele- 
ments of and whose links are the transitions obtained by the 
invasion process just described. Figure [T]represents a typical 
assembly graph. From the viewpoint of statistical mechanics, 
Q is the phase space of our system. 

The assembly process becomes a Markov chain if to every 
pair of communities E and E' we assign the transition proba- 
bility Pee' = 5 EE' + £,Qee', where Qee' is the fraction of 
the L + 1 different invasions of E that lead to E' {Qee' 
provided Q contains the link E E'). P is the transition ma- 
trix of a finite, aperiodic, Markov chain, so its states are either 
transient or recurrent. There can be one or several subsets of 
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FIG. 3: Number of recurrent states of the Markov chain as a func- 
tion of the resource saturation 7?, obtained using the algorithm of 
Ref. (27I1 . Inset: Mean number of species in the stationary state vs. 
R. 



recurrent states, the chain being ergodic in each of them |12C 
Every recurrent subset is a different endstate of the assembly 
process. The assembly will depend on the history only if there 
are at least two such recurrent subsets. Each recurrent subset 
has a stationary probability distribution on its states. 

It is worth noticing that the computation of P is exact. This 
means that we have a complete and exact description of the 
evolution of the assembly process on the phase space. In par- 
ticular this means that we can compute the evolution of any 
observable in an exact — albeit numeric — manner, without 
resorting to taking averages over realizations of the process. 
This is the most important difference of this model with re- 
spect to all assembly models considered so far 

Let us now describe the results. First of all, the process has 
a unique set of recurrent states for all < i? < 1700. This 
means that for this model we are able to prove that the end- 
state does not depend on the assembly history |5]. This result 
agrees with previous evidence found in other assembly models 
a si ll in as well as in experiments |3], where the same kind of 
assumptions about the invasion rate are made. For some val- 
ues of R this endstate contains a single community but typ- 
ically it contains more than one, often very many. Figure |2] 
illustrates one of these sets along with the transitions between 
its states. At any time the ecosystem is realized by one of the 
communities of the set. Invasions may induce transitions be- 
tween communities of this set but never lead outside it, hence 
the set as a whole is resistant to invasions. The frequency with 
which a community is visited defines a stationary probability 
on the set. Notice that only a few communities (rather similar 
to each other) are visited with a high probability, so the net re- 
sult is as if the community were "fluctuating". Figure[3]shows 
the number of communities in the endstate vs. R. One can 
see that as R is increased towards the onset of a new level ap- 
pearance the number of communities increases considerably, 
to drop down to just one once the new level has established. 

Inset in Fig.[3]represents the mean number of species in the 
endstate vs. R. The dependence is basically linear, except for 
some dips near the onset of a new level followed by a discon- 



4 



Q." 



l.Or 
0.9 

0.8- 
0.7- 
0.6 
0.5 



- % 


= 0.05 


- 5 


= 0.10 


..... ^ 


= 0.20 




= 0.30 



0.3 — 
iU> 0.2- 

^"o.i-/ 
o.oj/- 



100 200 



300- 



0.4, 



"0 



20 



40 




- 5 


= 0.05 


- ^ 


= 0.10 


-- ^ 


= 0.20 


4 


= 0.30 



O.OL 




100 



200 



300 



400 



500 



FIG. 4: Probability of invasion, Pi (t), vs. mean number of invasions, 
^f. Inset: probability of reconfiguration after invasion, Pait), vs. 
^t. Upper pannel: R — 430 (endstate of Fig. Lower pannel: 
R = 540 (endstate with a single community). 



tinuous jump once the new level is established. This behavior 
reflects a top-predator effect: Top predators control the pop- 
ulations of their preys, avoi ding their extinction by overcon- 
sumption of their resources ll28ll . As a matter of fact we have 
checked that when a new level appears it contains a single top 
predator and the number of species at lower levels rises. 

Other interesting observables are the probabilities of ac- 
cepting the invader, Pi{t), and of undergoing a reconfigura- 
tion after the invasion, Pa{t). Both are obtained as Pi/a{t) = 
^^(^'^, Pee'){P*') E0, where the inner sum runs over tran- 
sitions from E to E' E in which the invader is accepted, or 
in which there is a reconfiguration of the community, respec- 
tively. Figure |4] illustrates them for two typical cases: One 
with a complex endstate and another one in which the end- 
state is a single community. Resistance to invasion increases 
as the comminity assembles. Species richness can be com- 
puted similarly. It grows monotonically up to saturation in the 
steady state. Both features, increasing resistance to invasion 
and increasing species richness, are common results of stan- 
dard assembly models 

In summary, this minimalistic assembly model exhibits the 
same behavior as those reported in the literature, indicating 
that this behavior is very robust and probably shared also by 
real ecosystems. Our model, however, provides a complete 
and exact description of both, the set of microstates and the 
dynamical pathways of the assembly process. So we are not 
limited, as in standard assembly models, to compute averages 
over a small set of realizations of the process. To give a hint 
about what this means, we have calculated, for R = 300 (a 
case with an endstate made of a single community of three 
trophic levels and 50 species), that there are ~ 10^° differ- 



ent minimum-length pathways leading from the empty to the 
endstate. This number is nothing that a simulation can come 
close to. This model allows us to prove rigorously that its end- 
state does not depend on the assembly history. Whether this is 
a feature that real ecosystems exhibit will, of course, depend 
on how well they fulfill the assumptions about the invasion 
rate underlying this and other assembly models. But a caveat 
should be made: As species of each trophic level are indistin- 
guishable, uniqueness (and hence independence on history) 
refers only to the number of species at each level. 
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